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The description of extreme-mass-ratio binary systems in the inspiral phase is a challenging problem 
in gravitational wave physics with significant relevance for the space interferometer LISA. The main 
difficulty lies in the evaluation of the effects of the small body’s gravitational field on itself. To 
that end, an accurate computation of the perturbations produced by the small body with respect 
the background geometry of the large object, a massive black hole, is required. In this paper 
we present a new computational approach based on Finite Element Methods to solve the master 
equations describing perturbations of non-rotating black holes due to an orbiting point-like object. 

The numerical computations are carried out in the time domain by using evolution algorithms for 
wave-type equations. We show the accuracy of the method by comparing our calculations with 
previous results in the literature. Finally, we discuss the relevance of this method for achieving 
accurate descriptions of extreme-mass-ratio binaries. 
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I. INTRODUCTION 

Extreme-Mass-Ratio Binaries (EMRBs) in the inspi¬ 
ral stage of their evolution are considered to be a pri¬ 
mary source of gravitational radiation ELi to be de¬ 
tected by the proposed laser interferometric space an¬ 
tenna LISA 0,3 [a, Hi- They consist of a “small” object, 
such a main sequence star, a stellar mass black hole, or a 
neutron star, with mass m ranging from 1M 0 to 10 2 Mq, 
orbiting a massive black hole (MBH) with mass M rang¬ 
ing from 1O 3 M 0 (if we consider the case of intermedi¬ 
ate mass black holes in globular clusters) to 10 9 Mq (the 
case of big supermassive black holes sitting in the center 
of galaxies). This translates to EMRBs with mass ra¬ 
tios, = m/M , in the range 10 _1 — 1CU 9 . In order to 
exploit this type of systems through LISA, it is crucial 
to have a good theoretical understanding of their evolu¬ 
tion, good enough to produce accurate waveform tem¬ 
plates in support of data analysis efforts. Because there 
is no significant coupling between the strong curvature 
effects produced by the MBH and its companion, rela¬ 
tivistic perturbation theory is a well suited tool to study 
EMBRs. Clearly, the accuracy of this approximation de¬ 
pends on the smallness of the mass ratio [i. 

The challenge in modeling EMRBs is to compute the 
perturbations generated by the small body in the (back¬ 
ground) gravitational held of the MBH, and how these 
perturbations affect the motion of the small body itself. 
This problem has been known in literature as the radia¬ 
tion reaction problem. This is an old problem and several 
approaches to deal with it have been proposed (see the re¬ 
cent review by Poisson and the contributions to j|j). 
The most extended approach consists in modelling the 
small object by using a point-like description and then, 
to describe the radiation reaction effects on the dynam¬ 


ics as the action of a local self-force that is responsible 
for the deviations from geodesic motion. A consistent 
derivation of the equations of motion coming out from 
this set up was given by Mino, Sasaki and Tanaka 0, 
and later, adopting an axiomatic approach, by Quinn 
and Wald [l(| (see also |I3). However, these works only 
provide a formal prescription for the description of the 
orbital motion. For the practical calculations of the self¬ 
force some techniques have been proposed: the mode-sum, 
scheme and a regularization scheme based 

on zeta-function regularization techniques [lb] (see |]j| 
for a recent progress report). 

The computation of the self-force and waveforms, and 
any other physical relevant information related to the 
inspiral due to radiation reaction constitute the main 
challenge of this problem. One possible way is to resort 
to analytic techniques by adding extra approximations 
to the problem, similar to those from post-Newtonian 
methods. However, the results may not be applicable 
to situations of physical relevance involving highly spin¬ 
ning MBHs and very eccentric orbits. To make computa¬ 
tions without making further simplifications of the prob¬ 
lem, numerical techniques appear to be a necessary tool. 
It is important to distinguish between frequency-domain 
and time-domain calculations. The frequency domain ap¬ 
proach has been used for a long time; it provides accu¬ 
rate results for the computation of quasinormal modes 
and frequencies (3I3l- However, the frequency-domain 
approach has more difficulties when we are interested in 
computing the waves originated from highly eccentric or¬ 
bits since one has to sum over a large number of modes to 
obtain a good accuracy. In this sense, calculations in the 
time-domain can be more efficient for obtaining accurate 
waveforms for the physical situations of relevance. 

However, the time-domain numerical approach has to 
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face a challenge, which consists of dealing with the differ¬ 
ent physical scales (both spatial and temporal) present in 
the problem and that expand over several orders of mag¬ 
nitude. Specifically, one needs to handle not only large 
wavelength scales comparable to the massive black hole, 
but also to resolve the scales in the vicinity of the small 
object where radiation reaction effects play a crucial role. 
The conclusion is that we need to incorporate adaptive 
schemes in our numerical algorithms in order to provide 
the resolution that every region in the physical domain 
requires. Since the small object is going to be moving 
through the domain (unless we choose a very particular 
coordinate system), it is convenient to allow the adap¬ 
tive scheme to change in time to distribute properly the 
resolution. Our choice to deal with these issues is the Fi¬ 
nite Element Method (FEM), which is a numerical tech¬ 
nique where adaptivity can be implemented in a natural 
way. The FEM has other properties, which we will dis¬ 
cuss in this paper, that make it very suitable to be used 
for the description of EMRBs and also for other physical 
systems that are the subject of investigation in Numer¬ 
ical Relativity. In a recent work |2(j, we have already 
tested Adaptive Mesh Refinement techniques intrinsic to 
the FEM in a toy model consisting of a particle orbiting 
a black hole in the context of scalar gravity, where we 
have shown how an adaptive scheme can provide better 
accuracy than a non-adaptive scheme with an equivalent 
computational cost. 

In this paper we use the FEM to perform time-domain 
simulations of a point-like object orbiting in geodesics 
(no radiation-reaction) around a non-rotating MBH, and 
compute physically relevant quantities like energy and 
angular momentum emitted in gravitational waves and 
waveforms. This type of calculations constitute a good 
touchstone to evaluate the Finite Element (FE) tech¬ 
niques that we present in this paper, specially in rela¬ 
tion to use this type of computations for the evaluation 
of the self-force on the particle. Since the MBH is non¬ 
rotating MBH the problem can be reduced to solve the 
one-dimensional Partial Differential Equations (PDEs) 
of black hole perturbations theory. These equations, in 
the Regge-Wheeler gauge, reduce to a master equation 
from which the metric perturbations can be fully recov¬ 
ered. The master equation for axial modes is known as 
the Regge-Wheeler equation, and for polar modes as the 
Zerilli-Moncrief equation. In this paper, instead of using 
the Regge-Wheeler function, we use a modification origi¬ 
nally proposed by Cunningham, Price and Moncrief r211 
that puts the axial modes on an equal footing with po¬ 
lar modes, as described by the Zerilli-Moncrief function, 
in relation to computing energy and angular momentum 
luminosities, and waveforms. 

The plan of the paper is the following: In section ITTI we 
summarize the main results from (non-rotating) black- 
hole perturbation theory that we need in our compu¬ 
tations, including the explicit form of the source terms 
coming from the particle energy-momentum tensor. As 
far as we know, the expressions we present here for the 


sources associated with the axial modes, described by the 
Cunningham-Price-Moncrief master function, are new. 
We also perform an analysis of the discontinuities in the 
master equations due to the Dirac delta distributions 
that the source terms exhibit. In section EH we describe 
the numerical framework. We use a FE discretization for 
the spatial domain and a Finite Difference discretization 
in time. We start with the discretization of the domain, 
consisting of dividing the computational domain into dis¬ 
joint subdomains (the elements). Then, we describe the 
FE functional spaces, which are finite-dimensional func¬ 
tional spaces used to approximate locally (at each ele¬ 
ment) our solution. The next step is the derivation of 
the weak form of the master equations, which is an inte¬ 
gral form. It is important to remark that FE algorithms 
are derived from the integral form of the equations, in 
contrast with other numerical techniques where the dif¬ 
ferential form is used to obtain a discretization. From the 
weak form, we obtain the spatial discretization by impos¬ 
ing the vanishing of the residuals of our equations, which 
basically means to impose the vanishing of the compo¬ 
nents of the equations with respect to a basis of func¬ 
tions constructed from the FE functional spaces. This 
process leads to a coupled system of Ordinary Differen¬ 
tial Equations (ODEs) which has a close analogy with the 
equations governing the behaviour of a system of coupled 
oscillators. A very important point in the discretization 
process is the fact that, because the FE formulation is 
based on an integral form of the equations, we obtain 
automatically a discretization of the sources containing 
Dirac delta distributions and its first derivative without 
having to resort to sequences of functions approaching 
in some limit the Dirac delta, we just use the properties 
of these distributions at the analytic level in the weak 
form of the equations. To solve the resulting ODEs we 
introduce a collection of evolution algorithms to solve the 
equations in second-order form and which have parame¬ 
ters that allow us to control the appearance of spurious 
high-frequency modes, which are common in systems like 
the one we are studying having a very localized source. 
We finish this section by discussing the structure of the 
mesh, in particular how adaptivity is implemented and 
how we can change in time this structure as the par¬ 
ticle moves. In section rm we discuss the performance 
of the FE numerical code we have developed and com¬ 
pare results regarding the computations of energy and 
angular momentum radiated with previous works in the 
literature, showing in this way the accuracy that this 
method is able to achieve. We conclude in section El 
where we discuss the convenience of using the FEM for 
the simulations of EMRBs in the light of the results of 
this paper and describe possible ways to proceed in the 
future to make this goal a reality. Finally, we have in¬ 
cluded two appendices: In Appendix El we summarize 
the geodesic equations of motion for the particle, and 
in Appendix El we briefly describe the Gauss-Legendre 
quadrature method for evaluating numerically some of 
the integrals that appear in the FE discretization of our 
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equations. 

The conventions that we follow throughout this work 
are: Greek letters are used to denote spacetime indices; 
capital Latin letters are used for indices in the time-radial 
part of the metric; lower-case Latin indices are used for 
the spherical sector of the metric. We use physical units 
in which G = c = 1. 


II. SUMMARY OF PERTURBATION THEORY 
FOR NON-ROTATING BLACK HOLES 


Then, any covariant derivative on the spacetime can be 
written in terms of the covariant derivatives on M 2 and 
S 2 , plus some terms due to the warp factor r 2 , which can 
be written in terms of w A ■ 

Metric linear perturbations of a spherically-symmetric 
background can be decomposed in scalar, vector and ten¬ 
sor spherical harmonics |4(j,[4l]]. The scalar spherical har¬ 
monics Y tm are eigenfunctions of the covariant Laplacian 
on the sphere: 

7 ab Y?™ =-1(1 + l)Y (m . (5) 


Perturbation theory of black holes has a long his¬ 
tory. It goes back to the seminal work by Regge and 
Wheeler 23 and later by Vishveshwara, Zerilli and Mon- 
crief mMM for non-rotating black holes and to 
Teukolskv j25l hfij for rotating ones. At present, in the 
case of non-rotating black holes the metric perturbative 
scheme is completely de velop ed and well understood at 
the linear level (see [ 23 , l28l lidl IsflL IsH . Is3 liidj for re¬ 
views). In the particular case of perturbations induced 
by an orbiting point-like object, which is the situation we 
are interested in this paper, there are a number of works 
on it Here, we summarize the 

theory using a particular formulation that makes the dif¬ 
ferent expressions involved compact and self-explanatory. 

We start from the perturbative splitting of the met¬ 
ric into the background, the non-rotating Schwarzschild 
black hole metric g®°^, and the small deviations h a p\ 


A basis of vector spherical harmonics (defined for l > 1) 
is 

y 'im \r£m nlm b\r£m (n\ 

a = 1 -.a ) = e a *b ! W 

where the Y^ m, s have polar parity (they transform as 
(— 1 ) J , like the scalar harmonics, under parity transfor¬ 
mations, and are also called even-parity type) and the 
S f 0 rn ’s have axial parity (they transform as (—l ) i+1 un¬ 
der parity transformations, and are also called odd-parity 
type). A basis of tensor spherical harmonics (defined for 
l > 2 ) is 

y -tnn _ rvtnn _ A 1) -tr£m 

ab — Y Tab 5 ^ab — Y :ab ' ^ * Tab 5 \*/ 

nUm n£m fo \ 

D ab = D (a:b) ? \°) 


ga/3 = g a/3 + K0 ■ (1) 

Due to the spherical symmetry, the background mani¬ 
fold is the warped product M 2 x S 2 , where S 2 denotes 
the 2-sphere and Al 2 a two-dimensional Lorentzian man¬ 
ifold. The metric can then be written as the semidirect 
product of a Lorentzian metric on M , g^, and the unit 
curvature metric on S 2 , that we call 7 a b- 


where the Y*™ , Z£^ have polar parity and the S have 
axial parity. 

We then split the metric perturbations h a p into polar 
and axial perturbations, h a p = h^p + /i^, and these can 
be expanded in the basis of tensor harmonics as 

'O = £'$". C = (») 

£,m £,m 


&ct{3 


Sab 0 \ 

0 r 2 7 ab ) 


( 2 ) 


Hereafter, x A denotes a coordinate system on M 2 and 
x a a coordinate system on S’ 2 ; r = r( x A ) is a function 
on M 2 that coincides with the invariantly defined radial 
area coordinate. In Schwarzschild coordinates we have 


g AB dx A dx B = - fdt 2 + f 1 dr 2 , / = 1 - 2 M/r. (3) 

A vertical bar is used to denote the covariant derivative 
on M 2 and a semicolon to denote the one on S' 2 , thus 
we have g AB \c = 7 ab-.c = 0. We can also introduce the 
completely antisymmetric covariant unit tensors on M 2 
and on S 2 , e AB and e a b respectively, in such a way they 
satisfy: e AB \ C = f-ab-.c = 0, e AC e BC = - 6 f , and e ac e bc = 
—(5° . It is useful to introduce a vector held variable for 
the gradient of r: 


w A = r 1 r ]A . 


( 4 ) 


where 


/ 0 q[S e a m \ 

\ * 4 m s e s) 


( 10 ) 


7 p ,£m 


]~£m \rtm 
n AB 1 


i~£m \r£m 
n A Y a 


* 


\K 


£m \r£m 
Y ab 


(^£m g£m ^ 


( 11 ) 

and where we use asterisks to denote the symmetry of 
these tensors. All the perturbations, (scalar), h 1 ™ 
and qtj' 1 (vector), and K (m , G em , and q[ m (tensor), de¬ 
pend only on the coordinates of the 2-manifold M 2 . 

The components of the energy-momentum tensor of a 
point-like object are given by 


T a0 = m 


/ 


dr 

7 ^9 


u a uds 4 [x — z{t)\ 


( 12 ) 
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where m is the mass, r denotes proper time, z(t) is 
the trajectory of the object, g denotes the metric de¬ 
terminant, and <5 4 is the four-dimensional Dirac density 
(f d 4 Xy/^gS 4 (x) = 1). We choose to decompose in har¬ 
monics the contravariant components T“^ [the decompo¬ 
sition for the covariant ones follows immediately]. In this 
way, the polar components can be described in terms of 
the following quantities 

Qil] = 8 tt [c ii}T AB Y ern , (13) 

Js 2 

| a— r 2 (■ 

<i4 > 

Q} m = 8nr 2 fdnT ab Y^, ( 15 ) 

Js 2 

Q L = 327rr4 |^f| j s d P ^ 2 > ( 16 ) 

where the bar denotes complex conjugation. The axial 
components can be described in terms of the quantities 

P ^ = WTT)J s d ? TAa§ * m ' < 17 > 

p im = Jdp T ab S e S , (18) 

To simplify the equations we choose to work in the 
Regge-Wheeler gauge: 


however, in this paper we will use the master function 
introduced by Cunningham, Price and Moncrief [21 ], fol¬ 
lowing the definition used in 0E3: 

®cpm _ _^_ AB ( im J.m\ Oil 

(e + 2)(l-l) V fe l A r aQb J ■ ( 21) 

One reason for using this function is to have the formu¬ 
las for the energy and angular momentum radiated to 
coincide with the ones for the polar modes (see below). 
In this respect, the formulation for axial modes is on an 
equal footing with the one for polar modes. For the po¬ 
lar modes, we use the well-known Zerilli-Moncrief master 
function: 

K tm + j (h%w A w b - rw A K fjj*) , 

( 22 ) 

where A = (£+2){£— l)/2+3M/r. The key point in intro¬ 
ducing these quantities is the fact that the Einstein per¬ 
turbative equations can be decoupled for them, and the 
remaining perturbative variables can be recovered from 
them. As it is well know, the equations for 4'°™ (or 
T)™) and are wave-type equations of the form: 


2 r 


- £m 


t{t +1) 


-d 2 + dl - v RW/ZM 


(r) 




cpm/zm 

tm 


n ^.cpm/zm 


(23) 


where r* is the so-called tortoise coordinate (r* = r + 
2M\n(r/(2M) — 1)). The potential for the axial modes 
is the Regge-Wheeler potential 

Vr(r) = ^Ue+D-—) , (24) 

\ r ) 

and the one for polar modes is the Zerilli potential 


h tm = G e m = o, qt™ = o. ( 19 ) 

but a fully covariant and gauge-invariant approach can 
be found in pkO] . 

The perturbative equations can be decoupled by intro¬ 
ducing certain combinations of the metric perturbations, 
which are gauge-invariant. For the axial modes, it is 
common to use the Regge-Wheeler master function 

*Z = ~l™ A <iT, (20) 


E/» = 


/ 


r 2 A 2 


2 / -i , 3 M\ M 2 / M\ 

2A f ( 1 + \e H-I + 18—— ( Ae -\ -) 

r ) r A 


r J 
(25) 

where Ae = (£ + 2)(£ — l)/2. The source term for axial 
modes is given by 


tm 


^ r AB -ptm 

£{£ + 1 ) A|B ’ 


(26) 


and for polar modes by 


oZM 

°tm 


j^WAQtm ^ ( 


(1 + A^)A 


W°g AbQ£™\C - ^W A W B Qt B ~ £Qjm 


1 

rA 


Af(A f — 1) 


M M 2 ' 
3(2A^3)- + 2W 


g AbQ£™ 


r 


(27) 


When we restrict ourselves to the case of a point particle, by introducing the energy-momentum tensor d into 
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the previous expressions, we find that the source term for 
both polar and axial modes has the following (singular) 
structure: 

S(t,r) = G(t,r)6[r - r p (t)} + F(t,r)5'[r - r p (t)], (28) 

where 5 is the one-dimensional Dirac delta distribution 
and S' its first derivative. The function r p (t) describes the 
radial motion of the particle in terms of the coordinate 
time t. In the polar case, the explicit expressions for the 
functions F(t,r) and G(t,r ) can be found, for instance, 
in (32, (38i). They are given by 

GZ(t,r) = a e (r)Y im (t)+b e (r)Y^(t) 

+ c e (r)Y£(t) + d t (r)Z e ”(t), (29) 


F%(t,r)=A t {r)Y tm {t), (30) 


where the f-dependence in the right-hand side has to be 
understood as: ( t ) = (9 = 9 p = 7r/2 ,<p = ip p (t)). The 
different functions of r are given by: 


ae(r) 


8nm f 2 j 6 M A 
l + X e rA 2 [— p ~E~ p 



(31) 


h( r ) 


167 TlTl L p f 2 r 
1 —(— E p r 2 A 


8nm L p f 3 




1 + A e E p r 3 A ’ 


(32) 


(33) 


nent of the four-velocity, which can be substituted by the 
expression given in equation lT&2ll . 

We have computed the sources for the axial modes for 
the case in which these perturbations are described by 
the Cunningham-Price-Moncrief master function, since 
we are not aware of any reference in the literature con¬ 
taining them. The result is the following: 


d(( r ) = —327 xm 


( 1 — 2 )! L p p 
{£ + 2)\E p r 3 ’ 


(34) 


GZ M (t, r) = u e (r)S e ™(t) + v e (r)S% , (36) 


M r ) = 


8irm f 



(35) 


where E p is the particle energy per unit mass, L p is the 
orbital angular momentum, and u r is the radial compo- 


where 


FZ M (t,r) = B t {r)S^{t ), 


(37) 


Uji(r) = 327 Tm 


{£ — 2)! / 2 L p / 
(£ + 2)\r*El \ 



(38) 


v e (r) = 327rm 


(1—2)! f 2 L p 

(£+2)! r 3 El 


B e (r) = 327rm 


(e-2y. f 3 l p 

{(■ + 2)! r El 



(39) 


It is worth pointing out that for circular motion, there is 
only contribution from the vector harmonics. 


The singular structure of the source terms in the mas¬ 
ter equations, for both polar and axial modes, shown 
in equation OBI implies the existence of discontinuities 
in the master function at the particle’s location. Out- 
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side the location of the particle the solution is smooth 
(assuming it was initially). One can then divide the one¬ 
dimensional spatial domain into two different and disjoint 
regions (which will change in time as the particle moves): 
The region to the left of the particle (r < r p (t)), and the 
region to the right of the particle (r > r p (t)). The master 
functions, on that regions, satisfy an homogeneous wave¬ 
like equation [Equation (1531) without sources]. Then, we 
can think of the solution of the full equation as being 
composed of the solutions of the homogeneous equations 
to the left and to the right of the particle, and relations 
between them at the particle location. Obviously these 
relations consist of imposing the discontinuities that the 
singular source terms dictate. In other words, we can 
write the solution of the full equation as 

«f(t,r) = 'f-{t,r)0(r p {t)-r) + 'f + (t,r)d(r-r p (t )), (41) 

where 6(r) is the Heaviside step function, and \k + (t,r) 
and 'F~(t,r) are solutions of the homogeneous equation 
to the right and to the left of the particle respectively. 
Due to the existence of the particle we will have that 


^ {t,r p ) i- V + {t,r p ), 


(a r T )(t,r p ) ± (d r ^ + ){t,r ). 

(42) 


By introducing m into the full equation we can derive 
the equation that govern the jumps in the master func¬ 
tions. The result is: 


(. f 2 (t)-r 2 p (t )) [0 r *]-2 r p (t)d t [4-] 

- (r p (t) + [*] = g(t, r p (t )), (43) 


are obtained via numerical integration of the geodesics 
ODEs shown in Appendix and the third one, r p (t), 
can be found directly from the geodesic equations: 

■■ ( f \ = (M. _ WA 

j r 2 P (t) El \ r p (t) 2 ) 

+ mm (i - . (49) 

We finish this section by giving the expressions of the 
averaged energy at angular momentum luminosities at 
infinity (as obtained from the Isaacson’s averaged energy- 
momentum tensor for gravitational waves |4J,[44|), which 
also hold at the horizon, in terms of the axial and polar 
master functions: 

E 7^(|*£T + I*?“| 2 ) , (50) 

l>2,m ' 

l = ^ E irn fr^. ■ 

n i>2,m ^ 

( 51 ) 

Finally, the metric waveforms are given by 

< " 2 ”' < 52) 

where _ 2 Y im are the spherical harmonics of spin weight 

—2 (see, e.g. 1451 1. 


(/ 2 M + r 2 p (t)) [T] = E{t, r p (t )), (44) 

where /(f) = f(r p (t)), and [T] and [c/.’F] are the master 
function discontinuities at the particle location 


= lim T + (f,r) — 

lim Ik (t,r), 

(45) 

r-»r p (t) 

r ¥r p (t) 

lim i9 r 4' + (t,r) 

lim 9T”(t,r), 

(46) 

r ^ r P (.t) 

r-»r p (t) 


where T and Q are functions of t and the particle loca¬ 
tion r p (t) that one obtains after applying the following 
properties of the Dirac delta distribution 

A(t, r)S(r - r p (t)) = A(t, r p (t))S(r - r p (t )), (47) 

A(t,r)6'(r - r p (t)) = ~{d r A)(t, r p {t))S(r - r p (t)) 

+ A(t, r p (t))S'(r — r p (t )), (48) 

to the original source terms iPl). The equations for the 
discontinuities [Equations fSll and (EH) ] contain the par¬ 
ticle radial position r (t), and its first and second time 
derivatives (r p (f), r p (t)). The first two, r p (t) and r p (t), 


III. THE NUMERICAL FRAMEWORK 

In this section we introduce the basics on the numerical 
framework that we use to solve numerically the pertur¬ 
bative equations described above. This task involves a 
number of choices that determine the particular features 
of our numerical method. To be specific, we want to 
solve our equations in the time domain , that is, we want 
to develop an algorithm that evolves the initial data from 
an initial state to a final time where we are interested in 
knowing the solution. This implies a discretization of 
our equations both in space and time. We choose to dis¬ 
cretize in space by using a Galerkin-type FE procedure, 
and in time by using Finite Differences techniques. In 
what follows we describe the details of these ingredients 
of our numerical calculations. 


A. The Mathematical Formulation of the Problem 

The model PDE problem that we are interesting in 
solving has the following form: 

£[’!'] = {-d 2 +dl~V(x))^{t,x)~S{t,x) = 0, (53) 
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^(t 0 ,x) = i/j 0 (x) , (d t ^)(t 0 ,x) = tj) 0 (x) , (54) 

[{d t -d x )^](t,x H ) = 0, [(Sf +d x )'S’]{t,x I ) = 0, (55) 
where 


t £ [t 0 , tf] , and x £ El = [x H , Xj\ . (56) 

Expression (15311 presents the structure of the equations 
we need to solve, namely a wave equation in a poten¬ 
tial V and with a source S. Here, x corresponds to the 
tortoise coordinate jy in the master equations, and for 
the sake of simplicity of the notation we will use it in 
most of the rest of the paper. The form of S(t, x) is as¬ 
sumed to be known, as it is in our case. Equation Q 
represents the initial conditions for the time evolution, 
we need the initial value of 4/ and its time derivative to 
solve a second-order hyperbolic problem. The boundary 
conditions, given in are simple one-dimensional out¬ 
going conditions (also known as Sommerfeld boundary 
conditions) at both ends of the spatial interval where the 
equation is considered, specified in m - This boundary 
condition is an exact outgoing boundary condition only 
at infinity, provided the potential and the source have a 
fall-off of the type of the potentials and sources that we 
are considering in our physical problem, otherwise it is 
just an approximate outgoing boundary condition whose 
accuracy depends on how far Xj is located. A similar 
argument also holds for the boundary condition at x H: 
the accuracy of the boundary condition that we use de¬ 
pends on how far towards —oo we locate x H . There is 
a well-known difference between x —> oo and x —► — oo 
which is due to the asymptotic behaviour of the poten¬ 
tials. When we approach —> oo the potentials behave 
like const, x x~ 2 , whereas when we approach —> —oo 
the potentials behave like const, x exp{— x/(2M)}, there¬ 
fore our approximate boundary conditions should work 
much better at x H than at Xj since our equations re¬ 
semble the wave equation better around x H . One could 
also get better boundary conditions by using the meth¬ 
ods suggested in references j4|| |^3| , where higher-order 
derivative boundary conditions are proposed, or by using 
the methods proposed in 14811491 l5Ct l5ll 1521 1 where exact 
radiative boundary conditions are studied. 


B. The Finite Element Discretization 

In this section we describe, in a simplified way, the 
main ingredients of the FEM that are relevant for our 
calculations. Detailed exposions can be found in 5jt, 54 

ii. 

The way one discretizes in space in a FE framework 
can be summarized in the following three steps: (i) Do¬ 
main discretization. The division of the spatial domain 


fl into a collection of disjoint subdomains {fl fc } i=1 N , 
the elements, (ii) The FE functional space. At every el¬ 
ement, f2fc, we introduce a finite-dimensional functional 
space, Fk, that we use to expand our fields locally at f Ik . 
Typically, these functional spaces are made out of poly¬ 
nomials. (iii) Weak formulation of the equation. This 
consists in converting the differential form of the equa¬ 
tions into an integral form that involves the boundary 
conditions of the von Neumann type, (iv) Equation dis¬ 
cretization. In a Garlekin-type of FE formulation, the 
discretized equations are obtain throught the imposition 
of the vanishing of all the residuals , £a = f n dxn t £[\I/], 
the components of our equation with respect a basis of 
nodal functions {n^x)} built out of the spaces Fk (see 
below). 


1. Domain Discretization and the FE functional spaces 

Points (i) and (ii) have to be treated jointly, because 
the structure of the elements and the structure of the FE 
functional spaces are not completely independent. 

Since we are dealing with a one-dimensional problem, 
the first step, the subdivision of the domain, is quite 
simple, we just divide the interval [x H , Xj\ into subinter¬ 
vals (see Figure m Hi = [x H , aq), n 2 = [x lt x 2 ), ... , 
Q n = [x N _ x ,Xj\ . This equivalent to locate N + 1 nodes 
in an ordered way: x 0 = x H < aq < • ■ ■ < x N _ 1 < x N = 
Xj . We denote the size (length) of the element il k by 
dk = x k +1 - x k ■ 


X N -3 x N-2 X N—1 X I 


d~Y ^ 


djv -2 


^N— 1 


djy 


FIG. 1: One-dimensional Mesh. 


The second step is very important in the sense that the 
accuracy and convergence properties of the FE scheme 
depend on the choice of the FE functional space. In this 
paper, we restrict ourselves to linear elements. For suf¬ 
ficiently regular meshes, linear elements lead to second- 
order convergence in the L 2 norm. On the other hand, 
the functional spaces Fk for linear elements are two- 
dimensional and can be span by the following two func¬ 
tions (see Figure |2l : 


M i{x) 


if x £ (xi,x i+1 ), 
0 otherwise, 


(57) 


Ni(x) 


Xi+ f if x £ (xi,x i+ i), 

0 otherwise, 


(58) 


which are usually called linear interpolation functions. 
However, to build a FE approximation of the solution of 







our PDEs it is more convenient to use the nodal functions, 
n t (x): 


FIG. 2: Linear interpolation functions M^x) and N t (x). 


n H {x) = N h {x ), n^x) = M N _ x (x ). (59) 

and for i = 1 ,N — 1: 

{ M i _ 1 (x) if x e 

N^x) if x G (xi,Xi+i), (60) 

0 otherwise. 

Their are called nodal functions because of the following 
property (see Figure 0: 

n i {x J )=6 ij , (61) 

that is, they vanish at all nodes excepting at the one they 
are associated with, where they take the unity value. 




*<-i *< * 


FIG. 3: Nodal functions n.i{x). 


2. The weak form of the equations 

The next step is to formulate our problem in what is 
called the weak form of the equation, which is an integral 
form. To that end, let us consider an arbitrary test func¬ 
tion </> on [t 0 ,tf] x [x H , Xj] and multiply equation 15311 by 
it. Then, we integrate over [x H ,Xj\ and apply integra¬ 
tion by parts to the term with second spatial derivatives. 
This produces a boundary term with first spatial deriva¬ 
tives that can be converted into first time derivatives by 
using our boundary conditions m - The result is: 


1 


£[4>, '&] = -! dx + (d x <j)){d x ^) + V{x)<j>'i\ + 4>d t ^\ x + 4>d t \V\ x -/ dx<j>S = 0. (62) 

Jx H 1 H Jx H 

I- 


This is the weak form of our equation, which has the re¬ 
markable property of including the boundary conditions 
of the problem. In the case we were dealing with Dirich- 
let boundary conditions, usually called essential bound¬ 
ary conditions in the FEM language, they would have not 
been incorporated in the weak formulation of the equa¬ 
tion. Instead, they are used to eliminate unknowns by 
providing values for them. 


3. The FE discretization of the equation 

To find the FE discretization of our equations we need 
to introduce the FE approximation to the solution of our 
problem. This approximation consists in an expansion in 


the nodal functions 


N 

Vh(t,x) =^ifj i (t)n i {x ), (63) 

i=0 


where the time-dependent functions ip i (t) are going to be 
the unknowns of our problem. Then, the FE discretiza¬ 
tion of the equations of our problem consists in imposing 
the vanishing of the residuals: 


S i = £[n i ,^h] = 0, i = 0, , N. 


(64) 
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This leads to one equation for node 1 , or equivalently, we 
have TV + 1 equations for our TV + 1 unknowns {ipf} • 
By introducing the FE approximation JMl into the weak 
form of our equation iPl) and imposing (IPl) . we get a 
system of coupled ODEs for the unknowns. We can write 
it in a matrix form as follows: 

M-vt> + (G-4' + IK-'4' = F, (65) 

where M, G, and K are (TV + 1) x (TV + 1) matrices, and 
\F and F are vectors with TV + 1 components, and IF 
and \F are the time derivatives of SF. The names and 
meaning of these objects can be thought to be inspired 
in the analogy of the system of equations lIPll with the 
system of equations for a coupled system of oscillators. 
The matrix M is usually called the mass matrix and its 
components are: 


our system of ODEs (l65l) analytically. The expressions 
for the components of mass matrix are: 

| [^i -1 j + 2 (di -i + d-i) + dA+i, 

(70) 

for i, j, = 1,. .., TV — 1, and the components related to the 
boundaries are: M hh = d H /2> , M h i = d H /6 , Mjv_i / = 
d N _ i/6, and M// = d N _ 1 / 3. The components of the 
damping matrix are simply given by 

= h ifSj H + $i idj i , ( 71 ) 

The first term of the components of the stiffness matrix 
is given by 


1 

di -1 


u i -1 


J 


+ 


1 

di -1 



1 

d i -1 


S 


i+lj i 


(72) 


r x i 

Mjj = / dxn i (x)nj(x ), (66) 

Jx H 

and hence it is a symmetric and positive definite matrix. 
The matrix G is usually called the damping matrix and 
has components 

= n^x^n^xn) + n^x^n^xf) , (67) 

which is symmetric and only has contributions from the 
boundaries, which means that our system only dissipates 
(or absorbs, if they sign of these two terms would have 
been negative, corresponding to ingoing boundary con¬ 
ditions) energy 2 through the boundaries and this reflects 
the fact that we are using outgoing boundary conditions. 
The matrix K is usually called the stiffness matrix. Its 
components are given by 

r x i 

K ij = / dx [n'^xjn'jfx) +V(x)n i (x)rij(x)] , (68) 

Jx H 

therefore, it is also symmetric. Finally, IF = (if>i(t)) and 
F is usually called the force vector and its components 
are given by 


r x i 

Fi(t) = — / dxS(t,x)n i (x). 
Jx H 


(69) 


for i,j,= 1,..., TV — 1, and the components related to 
the boundaries are: K hh = 1 /d H , K#i = —l/d H , 
Kjv-i/ = — l/djv_ 1 , and K// = —1 /d N _ 1 . The second 
term in the components of the stiffness matrix involves 
the potential, and therefore it has to be computed numer¬ 
ically. To that end, we use Gauss-Legendre quadratures 
(see Appendix iBl). 

We can compute the components of the force vector F 
by using the form of the source term S , which is given 
in equation ( 1551 ) . Using the properties of the Dirac delta 
distribution, we find that the structure of S implies the 
following structure for the components of the force vector: 


F i(t) = {v{r p {t))[{d r F){t,r p (t))-G{t,r p {t))] 

+ v\r p {t))F{t,r p {t))} n^Xpit)) 

+ v2 (r p (t))F(t,r p (t))rii(x p (t)), (73) 

where 


v{r) 


dx(r) 

dr 


1 

7’ 


v'(r) 


dv(r) 

dr 


2 M 

r 2 p ’ 


(74) 


and x p (t) is just the radial motion in terms of the tortoise 
coordinate. This completes the FE discretization of our 
problem. 


C. Evolution Algorithms 


In our particular case, we can compute most of the 
components of the matrices and vectors that determine 


1 In the case of problems involving essential boundary conditions 
we would get as many equations as unknows remain after impos¬ 
ing these boundary conditions. 

2 We can define an energy for our system by: 

E['E 4>] = i + ^ t -K -4-) . 

This energy would be preserved by the evolution when G = F = 
0. 


The next step is to solve the system of ODEs given 
in equation (OH), which is coupled to the ODEs corre¬ 
sponding to the motion of the point-like object [equations 
mm in Appendix m The numerical algorithms we 
use derive from the average acceleration method, which is 
based on the assumption that over a small time interval 
any nodal acceleration can be considered to be a linear 
function of time. Then, for a time interval (f 0 , t 0 + At), 
we write 

*(*) = *(tc) (i - + Ho + At)± . (75) 
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Integrating in time this equation twice and evaluating at 
t = ti = t Q + At we get 


9(to) + \ ^(t D ) + #(ti) 


A t. 


(76) 


V(ti) = Vito) + if (to) At + i 2#(t 0 ) + #(ii) (At) 


(77) 

This algorithm that one derives from these expressions is 
conditionally stable. Newmark introduced a gener¬ 
alization of the equations m and £73) in the following 
way 


*(t 1 ) = 4r(t 0 )+ (i- 7 )#(t 0 )+ 7 *(ti) At, (78) 


&{h) = *(t, 

l 


V(t 0 )At 


+ 21 


(l-2/3)*(t 0 ) + 2f3*(t 1 ) (At) 2 , (79) 


where ( 7 , (3) are parameters that have to be chosen for 
accuracy and stability. The Newmark method is uncon¬ 
ditionally stable for the following range of the parameters 

(7./3): 





(80) 


For ( 7 ,/J) = (1/2,1/6) we recover the average accelera¬ 
tion method; for ( 7 ,/?) = ( 1 / 2 , 0 ) we obtain the central 
differences method (although it is not strictly explicit), 
which is conditionally stable; and for (/?, 7 ) = (1/4,1/2) 
we get the trapezoidal rule, which is unconditionally sta¬ 
ble and second-order accurate. In the Newmark scheme, 
equation El is left untouched. 

However, numerical damping to prevent the amplifi¬ 
cation of high-frequency modes cannot be introduced in 
the Newmark algorithm without degrading the order of 
accuracy to first-order. There are a number of numerical 


schemes that generalize the Newmark scheme in order 
to include maximal dissipation of high frequency modes 
and minimal of low frequency modes and at the same 
time maintaining second-order accuracy. In particular: 
the Hilber-a method [13], the Bossak-a method |59|, and 
the Generalized-a method jSQ]. We present here the last 
one, which includes, for certain values of the parame¬ 
ters, the other methods. The Generalized-a method can 
be seen as a generalization of Newmark’s algorithm in 
the sense that equations m are also assumed by this 
evolution scheme. The generalization takes place when 
we discrete the set of ODEs given in equation El - Let 
(\F ra , \&„, be the values of our unknowns and their 
time derivatives at a time t = t n . Then, the discretiza¬ 
tion of El used in the Generalized-a method is given 
by 


M-# n+1 _ am + G-¥„ +1 _ a/ + K-¥ n+1 _ a/ = F n+1 _ af , 

(81) 

where 

i +1 - a „ = (l-«Ji +1 +aA, (82) 


*n+i-a, = (l-a f )¥ n+1 +a f ¥ n , (83) 


*«+!-«, = (! -a / )* B+1 + a / ^ n , (84) 


F n+1 -a f ~ (1 ~ a f)F n +l + a fFm ( 85 ) 

where a^ and a m are constants. The Newmark method 
corresponds to (aj,a m ) = (0,0), the Hilber—a method 
to a m = 0, and the Bossak—a method to a^ = 0. In¬ 
troducing equations El and El into equation El and 
rearranging the different terms, we arrive at the following 
equation for 


[(1 - a JM + (1 - a f )^AtG+ (1 - a f )/3( At) 2 K] -* n+1 = (1 - a f )F n+1 + a f F n - 


•T' 


-G- 


*„ + (!-«,)(!-7)A t* r 


+ (! - a f )At 


* n + (--/3)At* r 


( 86 ) 


Then, the algorithm that we are going to use to solve 
these equations for our unknowns goes as follows: (i) 
We solve El for (ii) We compute 'F n+1 from 

El and, (hi) We compute , ® , rl+1 from (1791) . Except¬ 
ing in very special cases, the method that comes out of 
this algorithm is implicit. In general implicit schemes 
are computationally expensive, but since we are using 


one-dimensional linear elements the matrices that we are 
dealing with are symmetric tridiagonal and therefore, one 
can use fast routines to invert them (see, e.g. 0 ). 

The convergence and stability properties of these al¬ 
gorithms and their high-frequency damping capabilities 
can be analyzed by casting the time discretization of our 
ODEs in the form: U n+1 = A-U n + R n , and then to 
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TABLE I: Values of the coefficients (a m , 07 ,/3, 7 ) that char¬ 
acterize the different evolution algorithms, in order to achieve 
consistency, stability, and favorable high-frequency mode 
damping properties. 


Algorithm 

ry 

'~ Aj m 

a f 

0 

7 

Newmark 

0 

0 

(1 + Poo) 2 

3 —Poo 

2(l+Poo) 

Bossak-a 

Poo -1 

P 00+ 1 

0 

l(l-a™) 2 

1 -a 

2 *-*ra 

Hilber-a 

0 

Poo 

3(1 + «/) 2 

h + a f 

!+Poo 

Generalized-a 

2 Pno — l 

Poo 

7(1 — Qm + Q /) 2 

2 H” Q f 

P 00+ 1 

1+Poo 


analyze the truncation error when U n is substituted by 
the exact expression and the spectral properties of the 
so-called amplification matrix A (see, e.g. j54|). A quan¬ 
tity that plays an important role is the spectral radius 


Poo = lim p(A), (87) 

A t/T —►oo 

where p(A) is the spectral radius of the matrix A, At is 
the time step, and T is the vibration period of a generic 
mode of the system. For p^ = 1 there is no damping, and 
the lower p^ is, the bigger the damping gets. In tabled 
we show the values of the time integration parameters 
(a m , 07 , /3, 7 ) for optimal damping properties in terms of 
the spectral radius p^ (see 0 E 2 f° r details). 

The numerical method we use to integrate the ODEs 
for the motion of the particle, equations (IA4IA5II . is 
the Bulirsch-Stoer extrapolation method (' IfidLldi i as de¬ 
scribed by |65] and |b6j (see also [6lJ). 


D. Structure and Motion of the Mesh 


elements selected into two elements of equal size. For the 
case in which the particle is at a node, we divide into two 
a given number of elements, say p # , to the right and to 
the left of the particle. We repeat this a given number of 
times, say q, . When the particle is located in the interior 
of a given element, we do the same but with the elements 
to the right and to the left of the element where the 
particle is located. In addition, we bisect the element 
where the particle is located q . times. Then, the mesh 
is determined by the three parameters {N a ,p t , q m ). The 
total number of elements in the case where the particle is 
located at a node is: N T = N 0 + 2 p # q % , and in the case 
where it is located in the interior of an element is given 
by N T = N a + 2 p. ( 7 . + 2 — 1. We show an example of 

these constructions in Figure 0] 



rm Nodes from refinement q, = 3 


gjgjogt t 

Nodes from refinement q, = 2 


Nodes from refinement q m = 1 


rr 
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| | Nodes fror 

1 refinement q, = 3 


11 

t 



! | I 

t 

t t 


1 


Nodes fro 

n refinem 

snt q, = 2 




Nodes from refinement q, = 1 


Element containing the Particle 


In the numerical simulations we have carried out we 
have used different mesh structures, all motivated by the 
fact that the size of the particle is very small compared to 
the length scale of the black hole. The features that dis¬ 
tinguish these different mesh structures are the following: 
(i) Refinement. Whether the size of the elements changes 
along the mesh in order to increase the resolution around 
the particle, (ii) Particle’s location. Whether the particle 
is located at the position of a node or, in contrast, it is 
located in the interior of an element, (iii) Motion of the 
Mesh. Whether the mesh is changing in time (in such 
a way that the part of the mesh with more resolution 
always contains the particle) or it is static. 

In the case without refinement, we just divide the mesh 
into a given number of elements, say N a , with the same 
size: d i = d , for all i. In the case where the particle is 
located at a node, since the location of the boundaries is 
also given, at least one element must have a size different 
from d. The way in which we refine the mesh to increase 
the resolution around the particle is by dividing a certain 
number of elements in the proximity of the particle a 
certain number of times. Each time, we divide each of the 


FIG. 4: Examples showing the structure of the Mesh for 
(Pm,q») = (4,3): On the top we have the case of a Mesh 
where the particle is located at a node. On the bottom we 
have the case of a Mesh where the particle is always in the 
interior of an element. 


Since the particle is moving, it may be convenient to 
adapt the mesh so that the finest region is always around 
the particle. We do this by just translating the structure 
of the mesh but without modifying it, excepting for the 
elements containing boundary nodes, whose size we need 
to change so that after the movement, the mesh fits in our 
domain. In other words, the resulting mesh is the result 
of applying a translation to all nodes and the translation 
distance is just the distance the particle has moved. Af¬ 
ter the mesh has been moved, we need to find the new 
associated nodal functions by using the FE interpolation. 
This is a very simple way of implementing a moving mesh 
technique, but given that in this problem we know at ev¬ 
ery moment where the resolution is required we do not 
need anything more sophisticated, at least at this stage. 
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IV. RESULTS FROM THE SIMULATIONS 

We have designed a numerical code with the ingredi¬ 
ents described above. To check its performance we have 
carried out a number of different test. First of all, we have 
tested it with Gaussian profiles propagating in flat space 
(both from rest and with initial velocities), and also with 
Gaussian profiles scattering off the potential of axial and 
polar modes. To that end we used uniform meshes and 
the following evolution schemes: the average acceleration 
method, and the Bossak-a, Hilber-a, and Generalized-a 
methods with different values of p x . In all these tests 
we found stable and second-order convergent (in the L 2 
norm) evolutions. Deviations from second-order conver¬ 
gence were found to be of the order of 0.1% . 

When we introduce the point-like object we are intro¬ 
ducing source terms that contain a Dirac delta term and 
a term containing the first derivative of the Dirac delta 
distribution. These are singular terms, the first one in¬ 
duces a discontinuity in the first derivative of the solu¬ 
tion, whereas the second one induces a discontinuity in 
the solution itself (see the description of these discontinu¬ 
ities given above in Section |H|. This loss of smoothness 
of the solution with respect to the case without particle 
is quite significant. The presence of Dirac delta distribu¬ 
tions can still be handled by the FEM without losing the 
convergence properties of the algorithms, but the inclu¬ 
sion of source terms with the first derivative of the Dirac 
delta is too severe to maintain the accuracy and conver¬ 
gence properties of the numerical algorithms (see f° r 
a discussion of this issue). As a consequence, the conver¬ 
gence in general drops from second to first order. There 
is however a way of preserving second-order convergence, 
consisting of locating the particle at a node and, instead 
of solving the equation with the force term due to the 
particle, we solve the equation without the force term 
in the region to the left and in the region to the right 
of the particle location using boundary conditions at the 
particle location that enforce the magnitude of the dis¬ 
continuities of the solution that the particle source terms 
dictate. The way of computing these discontinuities is 
to use the equations that govern their behaviour [equa¬ 
tions m and 63) ]. However, this way of proceeding 
has some drawbacks depending on the way we imple¬ 
ment it. Either it requires to change the structure of 
the matrices in the FEM discretization of the equation, 
transforming the linear algebra problem and making it 
similar to the one that we would get if we were using 
high-order elements, or it changes the stability proper¬ 
ties of our time-evolution schemes from being uncondi¬ 
tionally stable (they are implicit schemes) to be subject 
to a Courant-Friedrichs-Levy stability condition. More¬ 
over, locating the particle at a node means to change the 
mesh structure at every step in the evolution (excepting 
for circular orbits), which means to use the FEM inter¬ 
polation every single time step. The conclusion one can 
extract from this discussion is that each of the different 
possible ways in which we can carry out the computations 


have some advantages and some disvantages. The perfor¬ 
mance of each of these possible computational schemes 
depends on the physical setup we want to simulate and 
therefore, one has analyze on a case by case basis which 
is the appropiate method to use. 

With regard to the choice of the time scheme, which in 
our framework is equivalent to the choice of the param¬ 
eters (/3,7, a m , a^), the numerical experiments we have 
performed tell us that the inclusion of the particle gener¬ 
ates high-frequency modes that corrupt the solution and 
therefore we have to choose p^ different from unity to 
damp those unphysical modes. In the case of the New- 
mark scheme, this means to loose second-order conver¬ 
gence, and therefore it is not the best choice. More¬ 
over, from our numerical experiments we observe that 
the Bossak-a scheme is the one that seems to work bet¬ 
ter in the sense that it is the scheme that damps the 
high-frequency modes in a more efficient way. The other 
schemes seems to require a lower p^ (higher damping) 
than the Bossak-a method for the same performance. We 
have also seen that the optimal value of p x to be used 
depends on the physical case we want to simulate, circu¬ 
lar orbits are the ones that require less damping whereas 
highly eccentric seems to require much more damping 
(for a comparable pericenter distance). It also depends 
on whether we move the mesh or not and on the resolu¬ 
tion we use. 

In order to further assess the capabilities of our compu¬ 
tational framework and its adequacy for the type of phys¬ 
ical computations we are interested in, we have compared 
our numerical simulation with different results in the lit¬ 
erature for different types of orbits (geodesics). The ini¬ 
tial data for the master functions is zero data, that is, 
'F(t 0 ,r) = = 0. This creates an initial burst 

of spurious radiation which, after suhcient time, leaves 
the computational domain. The global spatial resolution 
we have used in the simulation varies from Ax = 0.1 M 
to Ax = 0.02M, and the number of times that we re¬ 
fine around the particle goes from q m = 0 to q, = 10 
(see subsection mm . The physical observers or detec¬ 
tors of the gravitational radiation are located at a tor¬ 
toise coordinate in the range |r„| = 2000 — 2500M, and 
the boundaries are located at a distance in the range 
|rj = 4000 — 6000M. Regarding the time step, it is im¬ 
portant to remark that because our evolution algorithms 
are implicit and unconditionally stable we are not sub¬ 
ject to a Courant-Friedrich-Levy type condition on At 
(excepting in the case where we use the scheme in which 
we impose the discontinuities generated by the particle 
at a given node), which we have taken to be At = 0.1M. 

To begin with, we compare results for circular orbits 
with the frequency domain code by Poisson [63, (as 
quoted in 33), and with the time-domain calculations 
by Martel [33 (using a formulation based on the Regge- 
Wheeler gauge and solving for the master functions) and 
Barack and Lousto |39| (using a formulation based on the 
Lorentz gauge and solving directly for the metric pertur¬ 
bations). The circular orbits considered have a radius 
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pM with p = 7.9456 and our observer is located in the 
radiation zone at r* = 2000M. We compare results for 
the energy and angular momentum luminosities to infin¬ 
ity with the results of J67,68] and @ in Table[nJ and for 
the energy luminosities to infinity with the results of 39] 
in Table Itm 

We have also compared results for elliptic orbits with 
the frequency-domain calculations of Cutler et al. 69] 
and with the time-domain calculations of Martel [381 ]. 
We have considered two types of elliptic orbits with or¬ 
bital parameters given by (p,e) = (7.50478,0.188917) 
and (p,e) = (8.75455,0.764124). For these orbits we 
have computed the averaged energy and angular momen¬ 
tum luminosities. The average is taken over a certain 
number of radial periods and our observer is located at 
r* = 2500M. The results are shown in Table ITvl 

We also compare results for parabolic orbits with the 
time-domain calculations of Martel [3^|. This type of 
orbits have e = 1 and are only characterized by the peri- 
center distance, which is given by pM/2 with p > 8. As 
p approaches 8 the number of orbital periods (A</? p /27r) 
diverges and the motion shows the so-called zoom-whirl 
behaviour (see, e.g. @), meaning that for a radial pe¬ 
riod the particle orbits close to the MBH for a number 
of orbital periods producing a very characteristic signal 
(see Figure Q with a number of cycles that depends on 
how close p is to 8. They are therefore a good test bed 
for the numerical computations. In Table m we show 
the computations of the total energy and angular mo¬ 
mentum radiated to infinity (E°°,L°°) and into the hori¬ 
zon {E h , L h ) for parabolic orbits with p = 8.00001 and 
p = 8.001. For these computations, our observers are 
located at r* = —2500M and r* = 2500M. 

We finish this section by commenting on the waveforms 
obtained from our numerical computations. We have al¬ 
ready mentioned that one of the advantages of the time- 
domain approach is that it can provide reliable waveforms 
for a reasonable computational cost. We show that this is 
indeed the case by plotting the following components of 
the waveforms: 'F™ f° r circular orbits with p = 7.9456 , 
d'g™ for elliptic orbits with (e,p) = (0.764124,8.75455), 
and \F™ for parabolic orbits with p = 8.001 in Fig¬ 
ures EJjSj and [7| respectively. To achieve a high degree of 
smoothness in the waveforms, the damping of the spou- 
rious high-frequency modes in the evolution is crucial. In 
this sense, our simulations show that the evolution nu¬ 
merical algorithms proposed in this paper are suitable for 
the production of reliable waveforms. 


V. CONCLUSIONS AND DISCUSSION 

In this paper we have presented a new method for com¬ 
puting the gravitational radiation emitted by a point like 
object orbiting a non-rotating black hole. We have shown 
that the method is accurate by comparing it with previ¬ 
ous results in the literature obtaining an agreement with 
relative errors of the order of 1%, in many cases even of 



Time (M) 


FIG. 5: Component (£, m) = (2, 2) of the waveform corre¬ 
sponding to circular orbits (e = 0) with p = 7.9456 . 


the order of 0.1% or below. We also have shown that 
these numerical techniques provide sufficiently smooth 
waveforms, which is one of the goals of these calculation 
in relation with gravitational-wave data analysis efforts. 
These results together with the particular feature of the 
computational method presented suggest that it is a suit¬ 
able method to be use in self-force calculations for inspi- 
ralling EMRBs and in the posterior waveform calcula¬ 
tions at the next perturbative order. Our numerical cal¬ 
culations are based on the FEM and related techniques. 
The main features of the FEM that makes it suitable for 
the study of EMRBs, and perhaps also for problems that 
Numerical Relativity deals with, are the following: (i) 
Proper description of the Computational Domain. This 
is particularly relevant when we want to solve the per¬ 
turbative equations in a 2D or 3D setup (see [13), as it 
is the case if we want to consider a rotating Kerr black 
hole, which is the astrophysically relevant case. It would 
be also relevant for the study of black hole spacetimes in 
Numerical Relativity. In this scenario, the spacetime ge¬ 
ometry may involve holes (inner boundaries arising from 
black hole singularity excision) and we may wish to use 
a spherical-type outer boundary to allow gravitational 
radiation leave the domain smoothly. All these geomet¬ 
ric issues have usually caused a number of problems in 
Finite Differences techniques, but can be handled in a 
natural way using the FEM. In this respect, the FEM 
has already shown its capabilities in solving problems in 
other scientific areas that involve much more complicated 
domains than the ones we can face in General Relativ¬ 
ity. (ii) Imposition of Boundary Conditions. In close 
connection with the previous point, the underlaying phi- 
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TABLE II: Comparison of the computations of energy and angular momentum luminosities at infinity for circular orbits with 
p = 7.9456 with results obtained with the frequency-domain code by Poisson [?57[f68ll and the time-domain code by Martel j38| . 
They are calculated at r* = 2000M. The energy luminosities are expressed in units of (M/p) 2 and the angular momentum 
luminosities in units of M/p 2 . In square brackets we have included the absolute relative difference (rounded to the largest 
value). 
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TABLE III: Comparison of the computations of energy lumi¬ 
nosities [expressed in units of (M/p) 2 ] at infinity for circu¬ 
lar orbits with p = 7.9456 with results obtained in the time 
domain by Barack and Lousto {3f|. They are calculated at 
r* = 2000A7. In square brackets we have included the abso¬ 
lute relative difference (rounded to the largest value). 
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losophy in the FEM is that one should use the mesh that 
adapts best to the geometric characteristics of the prob¬ 
lem we want to solve. In particular, to the boundary 
conditions, since it is not equally simple and convenient 
to impose outgoing radiation conditions in a rectangular 
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FIG. 6 : Component (£, m) = (2,1) of the waveform corre¬ 
sponding to elliptic orbits with e = 0.764124 andp = 8.75455 . 


boundary than in a spherical one. This also has an im¬ 
pact when we perform the FEM discretization, since it 
is based on the weak form of our equations, which can 
have built in the the boundary conditions. In the case of 
problems in 2D or higher dimensions, if the boundary is 
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TABLE IV: Computations of the total average energy [in units of (M/p,) 2 ] and angular momentum luminosities [in units of 
M/p 2 ], < E°° > and < L°° >, for elliptic orbits. They are calculated at r* = 2500 M. We compare with the results obtained 
by Cutler et al. IFifl using a frequency-domain numerical code and by Martel [ 2 U using a time-domain numerical code. We 
consider two different types of elliptic orbits: Orbit A : ( p , e) = (7.50478, 0.188917) . Orbit B : (p, e) = (8.75455, 0.764124). 
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TABLE V: Computations of the total energy [in units of M/p 2 ] and angular momentum [in units of p~ 2 ] radiated, both to 
infinity () and into the horizon ( E H , L H ), in parabolic orbits (e = 1). They are calculated at r* = —2500 M and 
r* = 2500M. In square brackets we have included the absolute relative difference (rounded to the largest value) with respect 
the results obtained by Martel ]38) using a time-domain numerical code. 
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FIG. 7: Component (£,m) = (2,2) of the waveforms corre¬ 
sponding to zoom-whirl parabolic orbits (e = 1 ) with p = 
8 . 001 . 


natural (adapted to the problem), the implementation of 
the boundary conditions becomes trivial (see, e.g. poj). 
This has advantages even in ID problems, like the one 
we have studied in this paper, where the imposition of 
boundary conditions like Sommerfeld or von Neumann is 
simpler than in a Finite Differences framework. This pa¬ 
per illustrates this fact, (iii) Treatment of distributions. 
Many description of EMRBs treat the small body as a 
point-like object, which despite being somehow unnatu¬ 
ral in General Relativity, allows us to perform computa¬ 


tions in a consistent way. The consequence of having a 
point-like object is that the equations that we have to 
solve contain source terms where Dirac delta distribu¬ 
tions and its derivatives (up to second derivatives in the 
case we were solving the Teukolsky equation sourced by 
a point-like object) appear. To deal with this kind of 
distributions in a Finite Difference framework is not an 
easy task, and the different ways in which one can han¬ 
dle them involved not trivial a priori regularizations of 
the distributions. Instead, in the FEM, the fact that the 
discretization is based on the weak form of the equations, 
an integral formulation, is a key point. We can evaluate 
the integrals that involve Dirac distributions analytically 
by using the properties of the distributions, without the 
need of using any regularization of those distributions. 
A sample of this has been given in this paper, where we 
used the weak formulation of the problem to discretize a 
source term containing the Dirac delta distribution and 
its first derivative. Then, the type of discretization we 
would get is in this sense analogous to the one proposed 
by Price and Lousto [3(1 (U HU, where they also used an 
integral form of the equations to discretize them. There¬ 
fore, using the FEM provides an additional advantage for 
the study of EMRBs where the small object is treated as 
a point-like object, (iv) Adaptivity. This is a key ingre¬ 
dient for the simulations of EMRBs. The calculations 
presented in this paper do not necessarily require adap¬ 
tivity, but they are an excellent benchmark to test these 
techniques. However, for the case of rotating massive 
black hole adaptivity may be the only way of performing 
physically realistic simulations. The FEM is a natural 
choice to achieve the high level of adaptivity required, 
both in the construction of the mesh and later by using 
any of the robust techniques of mesh refinement available 
(see 0 and references therein). 

Apart from these specific reasons, there are other mo¬ 
tivations in favour of using the FEM. In this sense it is 
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important to mention that because the FEM is based on 
piecewise (polynomial) approximations, it lies on sound 
mathematical grounds (see, e.g. jFS tol). From the point 
of view of building numerical codes based on the FEM, 
it is important to emphasize the high degree of indepen¬ 
dence of the different ingredients of the FEM discretiza¬ 
tion process |53l (5J, [53, |53 , which makes it very suit¬ 
able for modular programming. In addition, the FEM 
has been widely used in many areas of scientific research 
and, as a consequence, a number of FEM packages and 
tools are available for scientific computation. 

There are a number of ways of extending this work 
in order to improve the computational framework in or¬ 
der to simulate EMRBs, and in particular to evalute the 
radiation-reaction effects. From the computational side 
we can introduce higher-order elements (by using FE 
functional spaces with higher-order polynomials), which 
will improve the accuracy of the computations. From 
the physical point of view, we can change the descrip¬ 
tion of the gravitational field, meaning the formulation 
of the perturbative scheme. In this sense, to compute 
the metric perturbations using the Lorentz gauge, as it 
has been recently proposed by Barack and Lousto f39f, 
appear to be a very convenient choice for a number of 
reasons (see [^3 for a detailed discussion). Among the 
advantages of this approach it is worth to mention the 
following ones: (i) Because one is working with pure met¬ 
ric perturbations the sources do not contain derivatives 
of the Dirac delta distribution, and hence the solution 
of the equations is continuous at the particle location, 
which will improve the accuracy of the computations, 
(ii) Moreover, in contrast with the computations in the 
Regge-Wheeler gauge, we do not need a metric pertur¬ 
bation reconstruction procedure (just algebraic computa¬ 
tions) to evaluate the self-force, (iii) The regularization 
procedures to obtain the self-force have only been given 
in the Lorentz gauge. It also has some disavantages: We 
need to solve a coupled system of equations instead of sin¬ 
gle wave-type equations, and there are constraints that 
need to be satisfied along the evolution. 

In the astrophysically motived EMRBs, the MBH is 
highly rotating and therefore it is desirable to be able to 
repeat these calculations by using the Kerr solution as 
the background spacetime. This is a more difficult prob¬ 
lem since it involves three-dimensional PDEs (or two- 
dimensional if we factor out the dependence in the polar 
angle). In this sense, it is important to mention that the 
FEM techniques that have been presented and used in 
this paper can be transfered to the higher-dimensional 
problem of computing Kerr perturbations. For the same 
reasons that have been pointed out before, a promising 
approach may be to solve for metric perturbations of the 
Kerr black hole in the Lorentz gauge. 
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APPENDIX A: MOTION OF THE POINT-LIKE 
OBJECT 


To complete the description of our physical prob¬ 
lem we have to introduce the equations of motion for 
the point-like object, which follows the geodesics of the 
Schwarzschild background black hole spacetime. Then, 
the four-velocity of the particle satisfies: 

u a V a u 0 = 0, u a = . (Al) 

dr 

The static and spherically symmetric character of the 
background imply the existence of first integrals of the 
motion (energy and angular momentum), and as it hap¬ 
pens in the Newtonian case, the motion takes place on a 
plane that, without lost of generality, we can take it to 
be the plane 9 = tt/2 ( u e = 0). Then, the equations of 
motion are equivalent to the following relations: 




• (A2) 


In order to obtain a well-behaved system of ODEs at the 
turning points of the radial coordinate (r = 0) we can 
use the following alternative quantity 


pM 

1 + e cos x ’ 


(A3) 


where e denotes the orbital eccentricity and p the semi- 
latus rectum , which can be used as alternative constants 
of motion to the pair ( E p , L p ). Then, the two equations 
we need to integrate to determine the position of the 
particle are: 


d\ {p — 2 — 2e cos x) (1 + e cos x) 2 Vp — 6 - 2e cos x 
dt ~ Mp 2 \J{p - 2) 2 - 4e 2 

(A4) 


dp (p — 2 — 2e cos x) (1 + e cos x) 2 
dt Mp 3 / 2 yj{p — 2) 2 — 4e 2 


(A5) 


APPENDIX B: GAUSS-LEGENDRE 
QUADRATURES 

The integrals in the second terms of are com- 
puted by using Gauss-Legendre quadratures (see, e.g. [Hj 
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and Ifill h Given a function F(x) we approximate its in- sociated with the zeros and given by 
tegral over the interval [a, 6] by 


dxF(x) 


b — a 
2 


N 

E W ? F 
1=1 


f a + b 



(B!) 

where Uj is the J-th zero of the Legendre polynomial 
P N (u) (it has exactly N zeros) and Wf 1 are weights as- 


W? 


2 

(i - u f)[p^K)f 


(B2) 


An N-point Gauss-Legendre quadrature integrates ex¬ 
actly polynomials of degree 2TV — 1. 
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